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fuU sky maps of the Stokes I, Q and U parameters in 5 frequency bands, span- 
ning 23 to 94 GHz, and supporting items, such as beam window functions and 
noise covariance matrices. The processing used to produce the current sky maps 
and supporting products represents a significant advancement over the first year 
analysis, and is described herein. Improvements to the pointing reconstruction, 
radiometer gain modehng, window function determination and radiometer spec- 
tral noise parametrization are presented. A detailed description of the updated 
data processing that produces maximum likelihood sky map estimates is pre- 
sented, along with the methods used to produce reduced resolution maps and 
corresponding noise covariance matrices. Finally two methods used to evaluate 
the noise of the full resolution sky maps are presented along with several rep- 
resentative year-to-year null tests, demonstrating that sky maps produced from 
data from different observational epochs are consistent. 

Subject headings: cosmology: cosmic microwave background — instrumentation: 
detectors — space vehicles: instruments 



INTRODUCTION 



The Wilkinson Microwave Anisotropy Probe {WMAP) is a Medium-class Explorer mis- 
sion designed to produce full sky maps of the cosmic microwave background (CMB) radiation 
in five frequency bands centered at 23, 33, 41, 61 and 94 GHz. WMAP was launched on 2001 
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the design of the observatory and the results of the 



first year's observations have been previously published. Analysis of the first three years 
of W MAP data has now been completed and the results are presented in companion pa- 
pers JPage et allbood : iHinshaw etaPbood : ISpereel et al1l2006h . 



One of the major design goals of WMAP was the careful control of systematic er- 
rors to allow the production of high quality maps of the microwave sky. Systematic errors 
in the first year results were analyzed in several papers accompanying the data release. 
A nalyses of the beam shapes, beam window functions and associated errors were presented 
(I20031). Details of the data processing and related errors were discussed in 



m 



Page et al 



Hinshaw et al.l (j2003al ). Radiometer performance and systematic errors were analyzed in 
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Jarosik et al.l (l2003br) and iHinshaw et al.l (l2003a| ) while errors related to sidelobe pickup 
were explored in iBarnes et al.l (120031 ) . 



This paper updates and extends the previous analyses through the use of three years 
of observational data and addresses additional issues, such as data processing specific to the 
production of polarization maps. The overall instrument performance and improvements to 
the instrument modeling are described in § O Included in this section are updates to the 
radiometer gain models and beam window function analysis. In § [3] changes and additions 
to the data processing procedures, including the technique used to produce the maximum 
likelihood sky maps for Stokes I, Q and U parameters, and the evaluation of the pixel-pixel 
inverse noise matrix are described. In § |4] tests on the sky maps, including year to year 
comparisons and evaluation of map noise levels, are discussed. 



2. INSTRUMENT PERFORMANCE AND MODELING 

2.1. Observatory Status and Observing Efficiency 

The WMAP spacecraft and instrument continued to perform normally throughout its 
second and third years of science data collection, spanning 2002 August 10 through 2004 
August 9. Two station-keeping maneuvers were performed during each of the second and 
third years of WMAP observations. Each of these maneuvers resulted in the loss of approx- 
imately 10 h of science data, approximately 0.5 h for the maneuver itself and the remainder 
for the recovery from the thermal disturbance to the observatory. The observatory also en- 
tered a safe-hold mode on 2003 10 Aug, believed to be the result of a cosmic ray event. This 
occurrence resulted in a loss of approximately 60 h of data, caused both by the time spent 
in safe hold mode itself, and the time required for observatory temperatures to re-stabilize 
after the associated thermal disturbance. These periods, plus a short (~ 3 h) data loss due 
to a data recorder overflow caused by ground station difficulties that occurred on 2003 May 
25, resulted in the loss of a total of 103 h of science data, resulting in an overall observing 
efficiency of 99.4 % for WMAP's second and third years of observations. A tabulation of the 
times of the data excluded from processing by the aforementioned events can be found in 



Limon et al.l (120061 ). 



During the second and third years of operation, 8 sudden jumps were observed in the 
outputs of the radiometers. These events, distributed among 5 of the 20 radiometers com- 
prising WMAP, are believed to be the result of sudden releases of thermally ind uced mechan- 



ical st resses that slightly alter the properties of some radiometer components (jJarosik et al. 



2003bl ). During the first year of operation, 21 such events were observed. It is believed that 
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the reduction in the number of events is the result of a more stable thermal environment 
during the second and third years of operation. Although the events are of short duration, 
the data processing used to produce the sky maps requires that a 2-4 h segment of data for 
the affected radiometer be excluded from the final processing to avoid artifacts in the sky 
maps. All 8 events taken together caused a loss of 0.015 % of the science data from the 
second and third year of observations. 



2.2. Instrument Thermal Environment 



Providing a stable thermal environment for the WMAP instrument is a key component 
in the production of a high quality data set with well defined systematic error limits. Of par- 
ticular importance are temperature variations synchronous with the 129 s spin period of the 
observatory. The methods used to design and predict the thermal env ironment of the instru- 
men t and how they affect the instrument performance are described in lBennett et al.l (j2003d ) 
and iJarosik et al.l (j2003al ). On orbit, 56 high-resolution platinum-resistance thermometers 
monitor the temperature of key instrument components. These data provide verification that 
the actual thermal stability meets design specifications and allow for estimates of spurious 
signal levels resulting from radiometer thermal fluctuations. Analysis of the flight thermal 
data indicates the expected annual temperature modulation arising from the eccentricity 
of WMAP^s orbit and the asymptotic gradual warming of the instrument components as 
the thermal blankets degrade from ultraviolet exposure. Values characterizing the annual 
temperature modulation and gradual warming of the instrument components are presented 
in Table [H The slow warming of the instrument components has increased the radiometer 
noise levels by < 0.1 % yr~^ based on analysis of the noise in the sky maps as presented 
in § 14.31 Projections indicate that the WMAP observatory should be able to operate for 
years (with rare solar eclipse avoidance maneuvers) . Plots of the instrument and observatory 
temperature histories are presented in iLimon et al.l (120061 ). 



Significant spin-synchronous temperature variations are detected in 10 of the 13 ther- 
mometers attached to the thermal reflector system (TRS) with peak-to-peak variations rang- 
i ng from 15 — 198 /iK rms. Similar v ariations were seen during the first year of observations 



(jPage et al.ll2003l : iLimon et al.ll2003l ) and were identified as the result of impulsive heating by 
solar radiation scattered off the rim of the solar array. Figure [T] displays the largest of these 
signals for each year of observations. The measured temperature profile as a function of Solar 
azimuth with respect to the observatory follows the same pattern for all 3 years of data. The 
similarity of the temperature profiles indicates that no significant degradation has occurred 
in the solar screening properties of the solar arrays or insulation blankets. The predicted 
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microwave emissivity of the reflector surface (0.0005) yields an estimated spin synchronous 
peak-to-peak spurious signal of ~ 0.2 /iK in the time ordered data (TOD). 

No spin-synchronous temperature variations were detected in the Focal Plane Assembly 
(FPA), which contains the cold radiometer components, the Receiver Box (RXB), which 
contains the warm radiometer components and phase switch drivers, the Analog Electronics 
Unit (AEU) containing the data collection electronics, or the Power Distribution Unit (PDU) 
which contains the power supplies for the HEMT amplifiers. The limits for such signals are 
presented in Table [H The bounds are limited by thermometer readout noise and are similar 
to the year-1 values. Limits on possible spurious signals arising from thermal fluctuations 
of the radiometers are obtained by combining the on-orbit temperature variations measured 
during the second and third years of WMAP observati ons with instrument thermal sus- 
ceptibility coef ficients based on prefl ight measurements ( jjarosik et al.ll2003bl ) and on-orbit 
measurements (IHinshaw et al.ll2003al ). Artifacts in the TOD resulting from thermal fluctua- 
tions of radiometer components are found to be less than 0.4 /xK rms for all 20 radiometers, 
consistent with the year-1 results. Given the low level of the expected spin synchronous sig- 
nal, no corrections for spin synchronous effects have been applied to the time ordered data 
in the 3-year analysis. 



2.3. Pointing Determination 

The boresight directions of the WMAP beams relative to the observatory are measured 
by relating radiometric measurements of Jupiter to the attitud e of the observatory as de- 



termined by two onboard star trackers (IHinshaw et al.l l2003al ). WMAP observes Jupiter 



~ 90 days per year in two ~ 45 day "Jupiter seasons" . To test pointing stability, boresight 
directions are computed separately for each Jupiter season, and differences between the in- 
dividual seasonal determinations are noted. During the first year of WMAP observations, 
the azimuthal beam positions found for the first two Jupiter seasons agreed to better than 
3", but the elevation positions differed by ^ 10". This small difference was consistent with 
expected error in the spacecraft quaternions and was treated as part of the error budget 
rather than being actively corrected. 

The second year of WMAP operation provided two additional Jupiter observing sea- 
sons, 2002 Nov 01 - 2002 Dec 26 and 2003 Mar 13 - 2003 May 08. With the addition of the 
third season, it was found that, while azimuthal positions were stable, the apparent eleva- 
tions of the beams relative to the spin axis of the satellite now differed by ^ 30" from the 
positions originally computed from the first Jupiter season. Jupiter observations from season 
4 confirmed the systematic change in elevation. By comparing the attitude measurements 
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Table 1. WMAP Year-2 and Year-3 Thermal Characteristics 



Component 


Mean 


AT 

avg 


Annual Mod 


Meas. SS. 


Req. SS. 




(K) 


(K) 


(K) 


{liK rms) 


(/xK rms) 


Power Distribution Unit 


297.4-300.0 


0.98 


0.97 


< 6.9 


< 10000 


Analog Electronics Unit 


300.3-307.4 


0.87 


0.80 


< 8.3 


< 10000 


Receiver Box 


287.1-289.2 


0.44 


0.54 


< 8.0 


< 500 


Focal Plane Assembly 


88.2-91.6 


0.19 


0.16 


< 5.7 


< 500 


Thermal Reflector Assembly 


67.6-87.1 


0.16 


0.16 


~ 200 


< 5000 



Note. — This Table summarizes the thermal environment of the WMAP observatory 
during its second and third years of observations. The ranges specified for the "Mean" are 
the 24 month averages of the hottest and coldest thermometer in each assembly. The av- 
erage increase in temperature year-3 relative to year-2 is given as ATavg. "Annual Mod" 
is the amplitude of the annual temperature modulation arising from the eccentricity of the 
WMAP orbit. The measured spin synchronous (Meas. SS.) values are obtained by binning 
measured temperatures in spin synchronous spacecraft coordinates relative to the Sun, af- 
ter spline removal of long term temperature drifts. The values presented are derived from 
one month of data and are the largest limits obtained for any thermometer within each as- 
sembly. No significant spin synchronous temperature variations are observed for the Power 
Distribution Unit, Analog Electronics Unit, Receiver Box or Focal Plane Assembly, so the 
limits presented (indicated by < ) are determined by the binned readout noise of the ther- 
mometers. Spin synchronous temperature variations are detected in the Thermal Reflector 
System and are at levels similar to those observed during the first year of operation. The 
last column contains the design requirement for the spin synchronous temperature stability 
of each assembly. The measured spin synchronous temperature fluctuations are all far below 
the maximum design values. 
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Fig. 1. — The temperature of the middle of the B side primary mirror vs. the Sun azimuth 
in the spacecraft coordinate system. The red data were taken during the first year of obser- 
vations, the black during the second and the blue during the third. The similarity between 
the three indicates that the screening of the cold optics from solar radiation is still effective 
during the second and third years of operation. 



- 8 - 



of the two star trackers and the radiometric Jupiter observations, the systematic variation 
was traced to a temperature dependent flexure occurring in the spacecraft structure, which 
altered the pointing of the star trackers relative to the microwave telescope's beams. 

Processing of the 3-year WMAP data assumes a linear dependence of apparent star 
tracker elevation with spacecraft temperature, and corrects the spacecraft quaternions ap- 
propriately. A new set of line-of-sight vectors describing the telescope beam positions for use 
with the updated quaternions is available with the data releas^. The resulting pointing cor- 
rections are small compared to the size of the beams (> 12') and therefore produce negligible 
changes to signals in the maps, but will slightly change the noise patterns since observations 
near pixel borders may have moved into adjacent pixels relative to the first year analysis. 
Based on the data from the six Jupiter seasons, residual pointing errors after application 
of the corrections are estimated to be < 10". A detailed description of the modeling of the 



pointing correction is contained in iLimon et al.l (120061 ). 



2.4. Radiometer Gain Model 

The year-1 analysis of the WMAP data employed a model that related small changes 
in the radiometer gains (typically ~ 1 %) to values o f radiometer houseke eping data based 



on a physical model of the radiometer performance (iJarosik et al.l l2003bl ). Each detector 
was modeled separately and used the measured values of the FPA temperature, TppA, the 
radio frequency (RF) bias powers on the detectors, V, and three parameters determined by 
fitting the model to the hourly gain measurements, obtained from the CMB dipole signal. 
The form of the model used in the year-1 analysis was 

V-V, 

G = a7f^ (^) 

-fFPA — -^0 

The parameter Vq characterizes compression (non-linear response) of the amplifiers and de- 
tectors. To characterizes the variation of the noise temperature of the input amplifiers with 
physical temperature, and a normalizes the overall gain. An example of th e gain model 



performance fit to the hourly dipole-based gain measurements can be found in iJarosik et al. 



( ]2003bh . 



When the gain models fit to the first year observation data were extended to the 3- year 
data significant deviations became evident. These deviations became apparent as the input 
housekeeping data to the model {V and Tfpa) spanned a larger range than was used to fit 



All data products in the current release are available through http://lambda.gsfc.nasa.gov 
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the model, due to the slow warming of the observatory. Re-fitting the original gain model 
parameters {Vq,Tq and a) greatly improved the agreement between the model and dipole 
derived gain measurements, but significant errors still were evident. An additional term 
was added to the model to account for the variation of the compression of the amplifiers 
and square law diode detectors with their physical temperature by allowing the parameter 
characterizing the nonlinearity to have a weak temperature dependence. The resulting form 
of the gain model is 

F- Vo-/3(rRXB- 290) 

(^ = (y 7f 7f , (2) 

FPA — 

where Trxb is the temperature (in Kelvin) of the RXB thermometer closest to the detector 
and P is an additional parameter to be fit. Figure [2] shows the dipole derived gain measure- 
ments for the V223 detector and the improved gain model for this detector, fit to the entire 3 
years of data. Also shown in this Figure is the original gain model, using model parameters 
fit to year-1 data only, but applied to the entire 3 years of data . The new model significantly 
improves the fit over the entire 3 year period. The first-year analysis implemented a small, 
time dependent weighting of the TOD using the denominator of equation ([1]) as a measure of 
the input referenced noise temperature of the HEMT amplifiers. The values of Tq and (3 are 
strongly coupled in the current gain model and are not independently determined with high 
accuracy. It is therefore not possible to use the denominator of equation ([2]) as a measure of 
the radiometer noise levels. The three year processing therefore assumes uniform radiometer 
noise within each year. We maintain the year-1 estimate of the absolute sky map calibration 
uncertainty of 0.5 %. 



2.4-1- The Gain Model's Effect on Measurement of the Quadrupole 

Comparison of the year-1 and 3-year gain models indicated that many of the year-1 gain 
models displayed small errors in the predicted gain roughly linear in time, with an error on 
the order of 0.3 % yr~^. For each DA^, difference maps were formed by subtracting maps 
processed with the original gain model from those processed with the improved gain model. 
These maps displayed a feature with a several micro-Kelvin quadrupolar component, an 
example of which is presented in Figure [31 This signal was found to have similar morphologies 
in many of the maps with its amplitude correlated to the difference between the mean slopes 
of original and improved gain models. No other multipoles show significant correlation with 
the mean slopes of the gain models. (Sinusoidal differences between the gain models with 



^ A Differencing Assembly (DA) is a pair of differential radiometers connected to the two linear polariza- 
tions of a set of telescope feed horns. 
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200 400 600 800 1000 
Time (days) 

Fig. 2. — Comparison of the hourly gain determinations (black) based on measurement of the 
CMB dipole to two different versions of the radiometer gain model. These data are for the 
V223 detector and the time range spans the three years of WMAP science data collection. 
The blue lines are the original gain model derived by fitting the initial 310 days of data. The 
light blue region, to the left of the vertical red line, indicates the time range used to fit the 
model. The dark blue region, to the right of the red vertical line, is this model as originally 
fit, applied to the remainder of the data. The orange line is the new form gain model fit to 
all three years of data. The updated gain model is a significantly better fit to the hourly gain 
measurements. Note that the difference between the models contains a component roughly 
linear in time. 
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Fig. 3. — Difference between the temperature slcy maps produced using two different gain 
models. Raw data from from the V2 differencing assembly for year-1 was processed both 
with the original (year-1) and the improved (3-year) gain models. The map projection shown 
is in ecliptic rather than Galactic coordinates. The observed quadrupolar feature arises from 
imperfect subtraction of the velocity induced dipole signal in maps processed with the year-1 
gain model. Similar features are observed in similarly constructed difference maps for many 
of the DAs. 
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Fig. 4. — Amplitude of the quadrupole difference between temperature maps produced using 
the original and improved gain models. The horizontal axis is the difference between the 
mean slopes of the gain models. Each black symbol represents one Differencing Assembly. 
The line is a fit to these points constrained to pass through the origin and has a slope of 
9.9 fiK/{% yr~^). The red symbols are estimates of the amplitude of residual quadrupolar 
artifacts in the temperature maps after application of the improved gain model. Note that 
these values are distributed around zero indicating that the phases of the estimated residuals 
are random and should not bias the measured sky quadrupole. 
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an annual period were also examined and showed no correlations with any multipole of the 
difference maps.) These small errors in the gain model resulted in errors in subtraction of the 
^ 3 mK dipole signal. When processed through the map-making pipeline this residual dipole 
signal yielded a systematic feature with a quadrupolar component common to many of the 
sky maps. Since the form of the residual gain error was common to many of the radiometers 
the resultant spurious quadrupolar features were aligned, resulting in a biased measuremen t 



of the quadrupole moments in both the auto and cross power spectra (IHinshaw et al.ll2003bl ). 
The alignment of this spurious feature was such that it partially canceled the true sky 
quadrupole signal. Correction of this systematic error accounts for a part of the small 
increase in the rep orted value of the qua drupole moment, £{£ + l)C2/27r, from 154 /iK^ in 
the year-1 release (IBennett et al.l l2003bl ) to 220 /iK^, the largest change being the r esult 



of the change of the estimator use to evaluate the quadrupole (IHinshaw et al.l |2006| ) . It 



should be noted that the current value is still far smaller than the mean ex pected value of 



1220/iK^ obtained from the best fit A-dominated cold dark matter models ( Hinshaw et al. 



20061). 



Estimates of residual quadrupolar contamination remaining in the maps processed with 
the improved gain model are obtained by multiplying a coefficient relating the amplitude 
of the spurious quadrupolar signal to the slope error in the gain model, 602/ Sg, by an 
estimate of the residual slope error in the improved gain model, A^f. The value of 6c2/Sg = 
9.9 /iK/(% yr^^) was obtained by a linear fit to the data contained in Figure H] constrained 
to pass through the origin. The residual slope error in the gain model, Ag, was estimated 
as the sum of two terms: a fit to the residuals of the gain model to the dipole based gain 
measurements, and a term to allow for a possible bias in the slope of the dipole based gain 
determinations. The value of this second term was obtained from simulations of the complete 
calibration algorithm including effects of far sidelobe pickup of Galactic and dipole signals. 
These simulations indicate that the hourly dipole based gain measurements have slope biases 
ranging from +0.03 to —0.04 % yr^^. The magnitudes of these biases were summed with the 
magnitudes of the residual errors between the dipole based gain determinations and the gain 
model to yield an estimate Ag. Values of A^f are all less than 0.23 % yr^^. The resulting 
estimates of the magnitude of possible spurious quadrupolar signals in the 3-year maps are 
plotted in red in Figure H] . An estimate of spurious contributions to C2 is obtained by 
calculating the quantity 

AC2 = 2 Ag"^ Cr"^ (3) 

where the average is over the Ka-W band DAs, and the nominal value for the quadrupole 
moment is taken as i{i + l)C2°™/27r = 236 /iK^. The resultant estimate of the uncertainty 
on the quadrupole moment arising from gain model errors, AC2, is 31 /xK^. 
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2.5. Time Ordered Data Power Spectra 

The outputs of the WMAP radiometers exhibit ex cess power at low freq uencies charac- 



terized by a "1/f knee frequency", /knee as described in iJarosik et al.l (l2003bl ). Production of 
sky maps involves application of a filter to the TOD, as described in §3.4[ These filters are 
derived from fits to the autocorrelation function of the TOD after removal of an estimated 
sky signal based on preliminary sky maps. The measured autocorrelation functions for each 
radiometer are parameterized as 

r AC, At = 1, 

NiAt) = l a + b\ogi\At\) + c[\ogi\At\)Y + d[\ogi\At\)f, 1< |At| < AW, (4) 
I |At| > AWx 

where At is the time lag between data points in units of samples, the parameters AC, a, 
b, c and d are determined by fitting to the autocorrelation data, and Atmax is the time lag 
at which the fit crosses zero, typically ~ 600 s. These fits were performed on a year-by- 
year basis to allow for gradual changes in the radiometer noise characteristics, even though 
the radiometer noise properties are very stable. An example of this procedure is presented 
in the top panel of Figure [5l which presents the measured autocorrelation data and the 
parameterized fit for the year-3 observations of the Wll radiometer. Tabulated values of the 
parameterized filter coefficients for all years and DAs are available with the data release. 



2.6. Transmission Imbalance Measurement 

The year-1 WMAP analysis included a measurement of a set of transmission imbalance 
factors, Xim, characterizing the different transmission of sky signals from the A side and B 
side optics into the radiometers. The response, d, of the WMAP radiometers to sky signals 
Ta and Tb may be written 



d OC (1 + Xim)TA - (1 



(5) 



where the factors parameterize a departure from ideal differenti al radiometer perfor - 
mance, specifically the level of response to common mode input signals (IJarosik et al.ll2003bl ) . 
For an ideal differential radiometer Xim = 0, and the radiometer exhibits no response to com- 
mon mode signal. 

The transmission imbalance factors are determined by fitting the raw TOD to templates 
composed of pure differential and pur e common mode signals originating from the CMB 
dipole anisotropy (IJarosik et al.ll2003bl ). The year-1 analysis was performed using 232 days 
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Fig. 5. — Noise and Filter Properties of the Wll Radiometer. The top panel displays the 
measured autocorrelation function of the Wll radiometer noise (black diamonds) and the 
parameterized fit to these data (red line) as described in § 12.51 The data point at At = 
with value 1 has been omitted for clarity. The remaining plots illustrate the steps used to 
form the N^^ filters used in the conjugate gradient map solution, described in § 13.4.21 The 
second panel displays the noise power spectral density obtained from a Fourier transform of 
the parameterized noise autocorrelation function, while the third panel shows the reciprocal 
of this function. The last panel presents the N^'^ filter function obtained via a Fourier 
transform of the reciprocal noise power spectral density, again with the data point at At = 
and value 1 omitted for clarity. 
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of data, and did not remove an estimate of Galactic emission and CMB anisotropy before 
performing the fits, since reliable sky maps were not available at the time the analysis 
was performed. The 3-year analysis improves on the previous analysis in two respects; 1) 
estimates of the Galactic and CMB anisotropy signals are removed from the TOD, leaving 
only the dipole signal, before fitting to the aforementioned templates, and 2) three years of 
TOD are used in the analysis. A fit to all three years of data was used to determine the 
values of the transmission imbalance factors, which are presented in Table [2l These values 
agree with those from the year-1 analysis to the degree expected, but are more accurate due 
to the improvements in the processing. 



2.7. Beam and Window Function Determination 

Knowledge of the optical beam shapes is of critical importance to the scientific interpre- 
tation of the data. The shape of the measured power spectrum, and through that the values 
of the cosmological parameters, depends on the beam window functions. In turn, the beam 
window functions are determined from the shapes of the beam profiles. Uncertainties in the 
beam profile at the -20 dBi level (30 to 40 dB down from the peak, near the noise level of 
the measurements) can influence the window function at the 1% level. This sensitivity has 
motivated an extensive investigation of the beams. 

The beam modeling we present here is based on six seasons of Jupiter data acquired 
over three years of observations. The modeling has been significantly enhanced and the 
approximations made for the first year analysis have been reexamined. With the new models, 
the beam profiles can be extrapolat ed below the nois e level of the maps, thereby obviating 



the cutoff radius, 6'rc, introduced in lPage et al.l (120031 ). It is found that on average the beam 



solid angles are systematically larger by 1% than the year-1 estimates, leading to a systematic 
decrease in the window functions at £ > 50, leading in turn to a systematic increase in the 
power spectrum. For example in the 200 < I < 800 range, the new combined V and W 
band window functions are 1.5% lower 1.5(t) than the same combination for year one, 
justifying the assumptions and treatment in the year-1 release. We retain nearly the same 
uncertainties on the window function that were given for year-1. 

The CO- and cross-polarization beam profiles for the A side are shown in Figure [61 
These update previous versions of the figures by including the cross-polar response and data 
further into the tail of the beams. The figures are based on pre-fiight data, though the 
general agreement between the projections and measurements of Jupiter is excellent. The 
predicted and measured cross-polar patterns are in general agreement. 
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Table 2. Input Transmission Imbalance Measurements of the WMAP Radiometers 



Radiometer 




Radiometer 




Kll 


0.0000 ± 0.0007 


K12 


0.0056 ± 0.0001 


Kail 


0.0035 ± 0.0002 


Kal2 


0.0014 ± 0.0002 


Qll 


0.0010 ±0.0003 


Q12 


0.0047 ±0.0006 


Q21 


0.0075 ±0.0012 


Q22 


0.0103 ±0.0007 


Vll 


0.0011 ±0.0003 


V12 


0.0027 ±0.0003 


V21 


0.0043 ± 0.0006 


V22 


0.0051 ±0.0021 


Wll 


0.0092 ± 0.0024 


W12 


0.0023 ±0.0010 


W21 


0.0103 ±0.0019 


W22 


0.0084 ± 0.0018 


W31 


0.0014 ±0.0015 


W32 


0.0045 ± 0.0010 


W41 


0.0208 ± 0.0034 


W42 


0.0219 ±0.0062 



Note. — Measurement of the fractional input transmission imbalance, Xim, obtained from 
the WMAP 3-year data. They were obtained via measurements of the radiometer responses 
to the common mode signal arising from the CMB dipole. All the values are small, neverthe- 
less corrections for this effect have been included in the map making algorithm. The values 
of the uncertainties are estimated from the variance of measurements from 3 single year 
analyses. 
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The cross polar response of the main beam can be parameterized as a combination of 
a rotation of the hnear polarization reference axis around the line of sight, and a coupling 
to the Stokes V component. The relative magnitudes of these two components depends on 
the phase of the cross-polar coupling. The two outputs from each feed horn are coupled to 
different radiometers, so WMAP is not directly sensitive to this phase. While neither of these 
components can produce an errant polarization signal from an unpolarized input signal, they 
do affect the polarization measurement. The rotation of the polarization reference axis from 
their design directions was measured to be less then 1?5 during ground testing. The effect 
of the coupling to the Stokes V component is a reduction in the sensitivity to the Stokes 
Q and U components, presuming no Stokes V signal is present on the sky. The largest 
cross polar response is in W-band (-22 dB), corresponding to, at most, a 0.6% reduction 
in polarization sensitivity relative to the temperature calibration. At the levels indicated, 
neither of t hese effects alter the polarization measurements significantly. More details can 



be found in iPage et al.l (120061 ) . 



The beams of the V2 DA are used to illustrate the salient aspects of the beam analysis. 
Figure [7] shows the beam profile for year-1 and the three years combined. It also shows the 
accumulated solid angle as a function of angle from the main beam axis. Figure M shows the 
window function derived from the profile. A number of features are evident. 1) The profile 
is stable and is the same between years. 2) There are contributions to the solid angle at the 
-20 dBi level (the forward gain is 55 dBi), near the noise floor of the measurement. 3) The 
new beam transforms are slightly different from the old ones. 

The beams are modeled using a code based on the DADRA (Diffract ion Analysis of a 



Dual Reflector Antenna) physical optics routines (IRahmat-Samii et al.lll995l ) . The deviations 
from an ideal beam shape can be explained as deformations of the reflectors. The shape of 
the primary is parametrized with a set of 122 Fourier modes on a rectangular grid and the 
shape of the secondary with 30 modes. Based on pre-fiight measurements of the cold optics, 
the surface correlation length is ~ 10 cm, which is just resolved by our basis functions. 
A conjugate-gradient least squares method is used to solve for the shape of the primary 
and secondary reflectors. The parameters of the fit are the amplitudes of the modes. The 
model simultaneously solves for all ten beam profiles and takes into account the measured 
passbands. For simplicity we average over polarizations. Though the input can be thought 
of as 10 Jupiter maps with 10^ net pixels, the fitting is done in the time ordered data to 
obviate complications associated with map pixelization. The actual fit is a multistep process. 
At first just a few modes on the primary are considered. After an approximate solution is 
found, more modes on the primary are added, and finally the modes on the secondary are 
added. The solution is annealed and the resolution of the model is adjusted as the fitting 
progresses. So far, we have only run the model on the A side. 
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Fig. 6. — Predicted (top) and measured (bottom) focal plane for the A side for the co- and 
cross-polar beams. The contours are spaced by 3 dB and the maximum value of the gain in 
dBi is given next to selected beams. The measurement was done in the GEMAC (Goddard 
ElectroMagnetic Anechoic Chamber) beam mapping facility at NASA/GSFC. For both the 
predictions and measurements, measurements at twelve frequencies across each passband are 
combined using the measured radiometer response. The difference between the predictions 
and measurements are due to reflector surface deformations with an rms of 0.02 cm. In 



flight, deformations are larger (jPage et al.ll2003l ) and require the modeling described in the 
text. In K and Ka bands, the cross-polar patterns are clearly evident at the predicted levels. 
For the other bands, the polarization isolation of the orthomode transducer dominates over 
the optical cross-polarization leading to a higher cross-polar gain. This beam orientation is 
for an observer sitting on WMAP observing the beams as projected on the sky. 
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Fig. 7. — Left: Profile of tlie A-side V2 band beam comparing year-1 (black), 3-year combined 
(red), the model of the beam (green), and the Hermite fit to the beam profile (blue). The 
noise level of the data is apparent at the -30 dB level. The year-1 cutoff radius was 1?8. The 
improvements to the window functions come from understanding the beam in the -25 dB to 
-30 dB region. Right: The accumulated solid angle as a function of radius for A-side V2. 
The red line in this panel corresponds to the model (green) in the left panel. The vertical 
line corresponds to the cutoff radius used in the year-1 analysis. For this release, the integral 
is extended to 2?5. 




Fig. 8. — The beam transform for the A side of the V2 DA, comparing first-year and 3- year 
results. The top panel shows the beam transform from both the first year (red) and 3-year 
(black) analysis — they are nearly indistinguishable at this level. The bottom panel shows 
the difference between the beam transforms. Note that the 3-year beam transform is lower 
than the first year transform ne ar i 200 by about 0.5 %. This corresponds to 1% in the 
window function as discussed in iHinshaw et al.l (120061 ). 



- 21 - 



The fit method works because of the wide range of frequencies, the high signal to noise 
ratio, the angular resolution of the measurement, and the large sampling of the focal plane. 
Related methods generally employ so me sort of sampling of t he phase of the wave front — 
e.g.. Out of Focus (OOF) holography (IB. Nikolic &: Hillsll2002l ) — which of course cannot be 
done with WMAP. 

The resid u als of the model are shown in Figure [91 These should be compared to Figure 3 
of iPage et al.l (120031 ) which shows an earlier version of the model for the A side. For the 
best fit, the reduced xl = 1-22 (with roughly 10^ pixels). In other words, the model fits to 
near the noise of the measurement with just 152 parameters. One measure of the accuracy 
of the model comes from the comparison between the modeled and the measured beam solid 
angle as shown in Table [31 Excluding W2, the rms deviation is 1.6%; the W2 solid angle is 
predicted to be 4.8% smaller than is measured. This gives us confidence that we understand 
the beams. However, as can be seen in the figure, there are small residual features that 
the model does not capture. These features occur mostly in the steep parts of the profiles 
corresponding to the deformations with the greatest length scales. Nevertheless, the model 
is successful at unveiling faint large scale components of the beam. For example, a -26 dB 
Q-band lobe located 1?2 below Q2 was predicted by the model before it was found in the 
data. The model's effectiveness lies in the fact that it uses measurements with a high signal- 
to-noise ratio at all frequencies to find the shape of the reflector and thereby predict the 
beam shape many beamwidths off the beam axis at the noise floor of the measurements. 

Figure [7| shows a comparison between the measured and modeled beam. The model 
captures the main features of the beam though there are some discrepancies at the -20 dB 
level. Most importantly, the model reveals that the beam profile continues to drop with 
increasing radius. This is simply a result of the ^ 10 cm correlation length of the surface 
deformation s. With the new model there is no longer a need to impose a cutoff radius as 
was done in [Page et al.l (120031 ). 



The window f unctions are comp uted from the symmetrized beam profiles following the 
Hermite method in [Page et al.l (|2003l ). although the method is modified to take advantage of 
the modeling. The range of spacecraft orientations corresponding to the observational data 
contained in each sky map pixel depends on the ecliptic latitude of the pixel. For pixels near 
the ecliptic poles the observations occur nearly uniformly for all rotations of the spacecraft 
around each beam line-of-sight, effectively symmetrizing the beam. At lower ecliptic latitude 
the range of rotations is reduced. The use of the symmetrized beam profile approximation is 
therefore very good near the ecliptic poles, but worsens at lower ecliptic latitude. The least 
symmetrization occurs in the ecliptic plane. For beams in the ecliptic plane window function 
errors arising from use of the symmetrized beam approximation are less than 1% for i < 600 



-22- 




Fig. 9. — Beams in the WMAP focal plane. The top panel shows the measured beams, 
the middle panel shows the beam model, and the bottom panel shows the residuals. In the 
top two panels, each beam is scaled to its maximum, set to dB (red), and then plotted 
logarithmically to a level of -40 dB (blue). For the bottom panel each beam's residual 
is shown linearly as 100(data — model)/beam peak. The scales are ±10% for K; ±5% for 
Ka,Ql,Q2; ±3% for V1,V2; and ±2.5% for W. 
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in V-band and W-band. In Q-band the errors are somewhat larger due to the larger beam 
ellipticities, with window functions error less than 1% for £ < 300. The uncertainties in the 
window function determinations arising from incomplete beam symmetrization are included 
in the final window function u ncertainties. More details regarding beam symmetrization can 



be found in iPage et al.l (120031 ) . 



Figure [7] shows one such fit for the symmetrized radial profile of the V2 A-side beam with 
a Hermite expansion of order 170. To remove some of the sensitivity of the window functions 
to the noise tail, hybrid beams are constructed. The hybrid method retains the high signal 
to noise observations of the main beam and replaces the measurement with the model in 
regions of low signal. For the A side, the hybrid beams are constructed by (a) scaling the 
model to measured peak height, (b) choosing a threshold level, -Bthresh (Table [3]), and (c) 
using the measurement above -Bthresh and the scaled model below it. The hybridization is 
done with the full two dimensional profile in the TOD before the symmetrization. Since a 
complete model of the B side does not yet exist, the A-side profile is translated and rotated 
to the B side focal plane, scaled to the peak of the B-side data, and interpolated to replace 
B side observations in the TOD when they are less than -Bthresh- For both A and B sides 
the Hermite fits of order n = 170 are made to the symmetrized beam profiles. Analytic and 
numerical models show that the assumption of symmetrization is sufficient except in Q band 



where a 4% correction in made for i > 500 (IHinshaw et al.ll2006l ) 



The windo w functions are av ailable with the data release. In general they follow those 



in Figure 4 of iPage et al.l (120031 ). At low 1, the windows for polarization are negligibly 
different than those for temperature. At £ ~ 500, the polarization windows can differ 
from the temperature windows by a few percent because the effective central frequencies 
for polarization are different from those for temperature. Thus, the effective beams for 
polarization are different from those for temperature. We do not take this effect into account 
yet, as it is negligible at the current levels of sensitivity. 

The uncertainty in the window function is ascertained by changing the criterion for 
merging the model and the measurements, changing the cutoff in the Hermite fit, comparing 
the Hermite and Legendre polynomial methods, and propagating the formal errors from the 
fit. The adopted uncertainties, along with a comparison to year-1, are shown in Figure [TUl 
The new uncertainties are very close to the year-1 uncertainties except at low I in K through Q 
bands where the new uncertainties are larger. A typical window function has an uncertainty 
of 2-3%. These uncertainties add in quadrature in the cosmological analysis. We anticipate 
that the uncertainty can be reduced a factor of two after the development of the B side beam 
model and with more beam maps of Jupiter. 



In addition to the work described above, two end-to-end consistency checks are per- 
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Fig. 10. — The uncertainties of the beam transform functions for all bands. The uncertainties 
of the window functions are double these. The cyan lines are the new uncertainties. The 
black line is the uncertainty for year-1. The inner and outer magenta lines are the formal la 
and 2(7 uncertainties for the Hermite fits. The thin red line indicates the fractional difference 
between the Hermite and Legendre polynomial beam transforms. 
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formed. 1) We check to ensure that within any frequency band, all beams yield the same 
value for the temperature of Jupiter. This tells us that the beam solid angles are consistently 
computed for both A and B sides. We postpo r ie a fu ll reanalysis of the Jupiter calibration 
and recommend using the values in iPage et al. J2003h . 2) We make separate maps of the sky 
for the A and B sides and compute their power spectra. We then show that the ratio of the 
power spectra is unity to within the noise limits. This shows that the A and B side window 
functions are consistent. 



2.8. Sidelobe Corrections 

All radio telescopes exhibit some degree of sensitivity to sources outside of the main 
beam. Such pickup is a particular concern for CMB experiments since the measured signals 
are small and foreground emission can be large. The total beam of a telescope is traditionally 
described as the sum of two components, the main beam and the sidelobes. Although the 
demarcation between the two components is somewhat arbitrary, the two components sum 
to the total beam response of the telescope. The WMAP raw time ordered data, and the 
maps reconstructed from this data, represent the sky signal convolved with the total beam 
response of the telescope. 

The exponential factor in the Hermite expansion of the symmetrized beam profiles (see 
§ 12. 7p and the order of the polynomial fit determine the maximum radius for which the 
Hermite expansion accurately models the beam pattern. Outside of this radius the Hermite 
expansion quickly decays. The radius at which the Hermite expansion decays, rn, provides a 
natural point at which to separate the main beam and sidelobe response: the window function 
derived from the Hermite expansion represents the beam for r < th and the sidelobes account 
for the remainder of the beams. Ideally for r < rn, the sidelobe map would contain residuals 
between the true beam pattern and that parameterized by the Hermite expansion. However, 
since the main beam is modeled to levels below the noise floor of the measurements we have 
no way to determine these residuals. We therefore zero the area of the sidelobe maps for 
r < th. The values of rn are presented in Table |H These sidelobe maps, when convolved 
with various input signal maps, are used to implement corrections for the sidelobe pickup 
in the TOD as described in § 13. 3[ These maps differ from those contained in the first year 
data release in two aspects: 1) The region within a radius of th around each line of sight 
direction has been set to zero, and 2) The region of the W-band maps derived from the 
preflight sidelobe mapping has been scaled down by ~ 5 dB, based on a re-evaluation of the 
calibration of those measurements. These sidelobe maps are available with the data release. 
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Table 3. Parameters of the A-side Beam Model 



Band 


-Bthresh (dBl) 


Below peak (dB) 


^meas /^mod 


K 


17 


-30 


0.999 


Ka 


17 


-32 


0.982 


Q 


18 


-33 


0.971 


V 


19 


-36 


0.999 


w 


20 


-38 


1.018 



Note. — "Below peak" is the forward gain minus -Bthresh- It is the level in dB below the 
peak at which the measurement is replaced by the model, l^meas/^mod is the measured beam 
solid angle ( from observations of Jupiter) divided by the modeled beam solid angle. 

Table 4. WMAP Beam and Recahbration Factors 



DA th Recahbration factor 



Kl 


6?1 


1.0151 


Kal 


4?6 


1.0047 


Ql 


3?9 


0.9971 


Q2 


3?9 


0.9975 


VI 


2?5 


1.0009 


V2 


2?5 


1.0011 


Wl 


1?7 


1.0043 


W2 


1?7 


0.9985 


W3 


1?7 


0.9985 


W4 


1?7 


1.0033 



Note. — This Table lists recahbration applied to the time ordered data to compensate 
for calibration biases introduced by ignoring effect of the beam sidelobes during the fitting 
of the calibration solution. Also presented are the maximum radii for which the Hermite 
expansions of the radial beam profile are considered accurate. 
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3. DATA PROCESSING 

The 3-year WMAP sky maps were processed as three independent data sets, each con- 
taining 1 year of observational data. The maps were later combined to produce various data 
products, including a three year combined map. The processing used to produce sky maps 
differs in several aspects from that used to produce the year-1 maps. The most significant 
change is that the 3-year maps, comprising Stokes I, Q and U components, are the maximum 
likelihood estimates of the sky map given the radiometer "1/f" noise characteristics. The 
year-1 maps were unbiased representations of the sky signal, but wer e sub-optimal in term s 



of noise since they were produced from a pre- whitened TOD archive (IHinshaw et al.ll2003af ). 
Other significant processing changes include applying corrections for sidelobe pickup to the 
TOD for all DAs, inclusion of an estimated polarization signal in the fitting of the hourly 
baseline solution, and elimination of a weighting factor based on a radiometer noise estimate 
originating from the radiometer gain model. ( see § 12.41 ) 

The current data processing flow is outlined in Figure [TTJ Note that in the Make Sky 
Maps block three sets of maps are produced for each year of data. The first of these are full 
sky maps produced at HEALPix r9. WMAP utilizes the HEALPijJfl pixelization scheme and 
designates map resolution with the notation r4, r5, r9 and rlO, corresponding to HEALPix 
Aside parameter values of 16, 32, 512 and 1024, and approximate pixel side dimensions of 
3?7, 1?8, 0?11 and 0?06 respectively. 

These full sky maps suffer from small processing artifacts associated with observations 
in which one of the telescope beams is in a region of strong Galactic emission while the 
other is in a low Galactic emission region. These artifacts arise from slight errors in the 
radiometer gain determination and pixelization effects, resulting in errant signals in some 
high Galactic latitudes regions used in CMB analyses. The year-1 map processing eliminated 
this problem by only updating the value in the pixel w ith the high Galactic e mission for such 



observations in the iterative map making procedure (IHinshaw et al.l l2003al ) . Adapting this 
asymmetric masking to the conjugate gradient processing used to produce the current maps 
is not straightforward. Instead, a separate set of r9 maps, termed spm maps, is produced 
using a symmetric processing mask. These spm maps omit all observations for which either 
beam falls into a high Galactic emission region. The processing mask used to identify the 
high emission regions excludes 5.7% of the pixels, and is available with the data release. 
These spm maps do not suffer the aforementioned problem, but contain no data in the high 
Galactic emission regions. Data from the full sky maps are used to fill the unobserved regions 
of the spm maps. Details of this process are presented in § I3.4.2[ 



see http: / /healpix.jpl. nasa.gov 
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The last set of maps listed in the Make Sky Maps block of Figure [TT] is produced at 
rlO and contains only spm intensity maps. These maps are used to evaluate the high-£ 
temperature power spectra and are produced at higher resolution to minimize pixelization 
effects. To speed computation these maps were produced using only intensity information, 
so no corresponding high resolution polarization maps are available. 

The implementation of the entire data processing pipeline has been exhaustively verified 
through numerous simulations of TOD with the non-idealities described in § 13.4.11 both 
with and without inclusion of instrument noise. In the case excluding instrument noise the 
pipeline has been shown to reproduce the simulated input maps to sub-nanoKelvin accuracy. 
Simulations with instrument noise have been shown to produce unbiased e stimates of the 



i nput s imulated sky maps. The subsequent discussion follows the notation of iHinshaw et al. 



( l2003al ) and updates the description of the data processing therein. 



3.1. SMOC processing and Archive Generation 

Processing of the data in the Science and Mission Operations Center (SMOC) and 
Archive Generation steps (see Figure [TTj) is virtually unchanged from that used in year- 
1. The only substantive change is that corrections for the thermally induced star tracker 
position errors ( § 12.31) are applied to the star tracker data in calculation of the pointing 
quaternions included in the raw TOD archive. 



3.2. Calibration and Baseline Fitting 

The calibration and baseline fitting procedures, used to calibrate the radio metric data 



based on the CMB dipole signal, are similar to those used in the year-1 processing (IHinshaw et al 



2003al ). This is an iterative process in which both calibration data and approximate sky map 
solutions are generated. Hourly estimates of the radiometer gain are made by fitting each 1 
hour segment of the TOD, corrected for estimated sky signals (excluding the CMB dipole), 
to the sum of a template dipole signal and a baseline term that accounts for very slow drifts 
in the radiometer outputs. The template dipole signal is derived from WMAP velocity and 
attitude information and the measured barycenter dipole from the first-year WMAP results. 
Improved sky maps are then generated using the updated baseline and hourly gain solu- 
tions, which in turn are used to produce improved baseline and gain solutions. The 3-year 
processing is an improvement over the 1-year processing in that it corrects the TOD for 
estimated sky signals based on the Stokes I, Q and U components before fitting for the gain 
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Fig. 11. — Schematic overview of the 3- yea r WMAP sky map p rocessing pipehne. Sub- 
stantive changes from the year-1 processing (iHinshaw et aLll2003al . Figure 1) are indicated 
in boldface. 
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and baseline, whereas in the year-1 processing only the Stokes I signal was removed from 
the TOD before the baseline and gain solution were fit. 

The WMAP scan strategy ensures that the Stokes I sky signal components average very 
nearly to zero over a 1 hour precession period. However, the orientation of the polarization 
axis of the feed horns and scan pattern can transform certain sky polarization signals into 
TOD signals with very long periods, extending to many hours or days. The baseline fitting 
routine forces the mean of the signal corrected TOD to zero on time scales of one hour and 
longer. Removing the polarization signal before performing the baseline fitting ensures that 
the polarization signals contained in the TOD remain unbiased. 

After the hourly gain and baseline solutions are converged they are used to fit the 
parameters of the improved gain model described in § 12.41 



3.3. Calibrated TOD Archive Production 

The production of a final calibrated TOD archive involves two major steps: first the 
gain and baseline are applied to the uncalibrated TOD using the parameters determined in 
the calibration and baseline fitting procedure. Two corrections related to sidelobe response 
are then applied to this initial calibrated archive. 

The first sidelobe correction simply removes an estimate of the signal arising from 
sidelobe pickup. This signal is calculated by convolving the Stokes I sidelobe response 
maps (j^ 12.81) for e ach D A with preliminary Stokes I sky maps using the technique of 



Wandelt fc Gorskil (120011 ). The input sky maps contain three components: 1) An estimated 
CMB + Galactic signal obtained from a preliminary set of maps produced without sidelobe 
corrections, 2) a fixed barycenter CMB dipole signal, and 3) an annually modulated CMB 
dipole signal from Earth's motion relative to the solar system barycenter. These corrections 
apply to the Stokes I signal only and are therefore applied equally to the data from both 
radiometers comprising each DA. This correction is applied to th e TOD on a sampl e-by- 



sample basis. The effects of polarized sidelobe pickup are small (IBarnes et al.l l2003l ) and 
are not corrected in this data release. 

The second sidelobe related correction compensates for a small calibration error intro- 
duced in the initial calibration process by multiplying the sidelobe corrected TOD archive 
by an overall scale factor. The uncalibrated TOD archive used to fit the gain and baseline 
solution contains signal arising from both the main beam and the sidelobe response of the 
telescope optics to the true sky signal. However, the template dipole signal used to fit the 
gain is based on ideal pencil beams in the boresight directions of each telescope beam. This 
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approximation can lead to biased gain measurements since it ignores the component of the 
TOD arising from sidelobe pickup of the sky signaL This effect has been studied through 
simulations in which simulated TOD, based on the sidelobe response maps, is input to the 
hourly gain fitting algorithm to determine the level of bias originating from the sidelobe 
signals. The bias is found to be relatively small and varies only slightly over the course of a 
year as the precession axis of the WMAP scan pattern sweeps across the sky. Table H] lists 
the values of these recalibration factors. A detailed description of this correction is presented 
in Appendix IM 



3.4. Maximum Likelihood Map Solution 



As described in lHinshaw et al.l (j2003af ). the WMAP time ordered data may be written 

as 

d = Mt + n, (6) 

where M is the mapping matrix that transforms a sky map, t, into discrete samples, d, and 
n is the noise added to each sample by the radiometer. The map making problem consists of 
generating an estimated sky map given d and the statistical properties of the noise. Taking 
the noise description as 

(n) = 0, (7) 
(nn^) = N, (8) 

it is well known that the maximum likelihood estimate of the sky map, t, is 

t = (M^N-^M)"^ ■ (M^N^M). (9) 

In the case of a single differential radiometer and sky maps comprising only Stokes I, d is a 
data vector comprised of Nt time ordered data elements, t is a sky map of Np pixels, M is 
an Nt X Np matrix and N"-*- is an A^^^ x Nt element matrix. The matrix M is sparse, each 
row corresponding to one observation, and each column corresponding to a map pixel. For 
an ideal differential radiometer each row contains a value of +1 in the column corresponding 
to the pointing of the A-side beam for a given observation, and a value of —1 in the column 
corresponding to the B-side beam. Given these matrices it is possible to produce maps 
from the WMAP data set using equation The second term on the right hand side of 
equation may be evaluated directly, while multiplication by the first term, (M'^N~-'-M)~^, 
may be performed using a conjugate gradient iterative technique. It is straightforward to 
generalize this technique to process polarization maps. 
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As described in lBennett et al.l ( l2003d ). the signal from each telescope beam (A-side and 
B-side) is separated into two orthogonal linear polarizations by orthomode transducers at the 
base of the feed horns. Each pair of feed horns is associated with two differential radiometers 
comprising the DA. One linear polarization from each feed horn is fed into one differential 
radiometer, while the orthogonal set of linear polarizations are fed to the second radiometer 
comprising the D A. Details of the alignment of the polarization axis and beam boresights 
were presented in iPage et al.l (120031 ) . 



The outputs of radiometers "1" and "2", di and d2, for each observation correspond to 
sky signals 

di = i(PA) + q(PA) cos 27a + u(pa) sin 27a - i(pB) - qbe) cos 27b - u(pb) sin 273 (10) 



and 

d2 = Kpa) - q(PA) COS 27a - u(pa) sin 27a - i{pB) + q(PB) cos 27b + u(pb) sin 27b. (11) 

Here i(pA), q(pA) and u(pa) are the Stokes I, Q and U sky signals at sky position pa- 
The polarization angle of radion aeter "1" with r e spect t o the sky reference direction for sky 
position Pa is 7a as described in iHinshaw et al.l (l2003al ). Variables with subscript B are the 
corresponding values for the sky position observed by the B-side beam. 

The mapping matrix M is generalized to include polarization processing by expanding 
it to 2A^t X 3A^p matrix where the number of rows has been doubled since there are two data 
values from each observation, and the number of columns has been increased to accommodate 
3 maps, corresponding to the three Stokes parameters. Each row of the polarization mapping 
matrix has 6 non-zero elements, with ±1 in columns corresponding to the observed pixels 
(A-side and B-side for each radiometer ) in the i map and ± cos 27a, ± sin 27a, ± cos 27b and 
± sin 27b, in appropriate locations in the columns corresponding to the q and u maps. Each 
observation is associated with 12 non-zero values of the mapping matrix that are distributed 
in two rows. The noise matrix is similarly expanded to account for the signals from both 
radiometers. The noise of the two radiometers is uncorrelated, so it may be described by 
the relations, 



(m) = (na) = 0, 
(nin2) = 0, 
(ninf ) = Ni, 
(nsn^) = N2. 



(12) 
(13) 
(14) 
(15) 



The full noise covariance matrix is simply a block diagonal combination of Ni and N2. 
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The noise from different radiometers was verified as uncorrelated (eq.[T3]) by evaluation 
of the cross correlation coefficient 

C„ = , <°-°-> (16) 

a/ (nini)(nin2) 

where rii and n2 were obtained from the TOD by subtracting an estimated signal based 
on a combined 3 year final sky map. For all DAs, except Kl, the measured value of the 
magnitude of C12 (at zero time lag) was —0.001 < C12 < 0. Similar results were obtained 
from simulations of sky signal plus noise which were constructed to have zero noise cross 
correlations. The small values of anti-correlation observed arise from noise in the common 
sky map used to remove the sky signal from the TOD of the two radiometers. For K-band, 
incomplete removal of the intense Galactic signal due to small gain errors (~ 0.1%) lead to 
a measured C12 value of ~ 0.007, similar to values obtained from the simulation which also 
included gain uncertainties. 



3.4-1- Radiometer N on- Idealities 

If the WMAP radiometers were ideal, solutions based on the mapping matrix described 
in the previous section would reproduce the maximum likelihood estimate of the sky signal. 
There are two known instrumental non-idealities that affect the maps and are treated in the 
analysis: bandpass mismatch and input transmission imbalance. 



Bandpass Mismatch 



The construction of the mapping matrix M presented in §3.41 implicitly assumed that the 
microwave frequency response of the two radiometers comprising eac h DA were identical . 
The WMAP radiometers have slightly differing frequency responses (jJarosik et al.l l2003bl ) 
which can cause signals with spectra different from that of the ca libration source (the CMB 
dipole) to be aliased into polarization maps (IBarnes et al.l 120031 ). This spurious signal has 
the characteristic that it appears as a difference between the response of the two radiometers 
comprising a DA, di — d2, but is not modulated with the polarization angles 7a and 7b 
as a true polarization signal would be, since it only depends on the intensity and spectral 
index of the source region. This effect can be treated in the map making procedure by 
considering the signal from the radiometers to originate from 4 source maps, the original 
3 Stokes parameters plus a spurious map, s. Including this extra term, the TOD is then 
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described by the relations 

di = i(pA)+q(PA) cos27a+u(pa) sin27A+s(pA)-i(PB)-q(PB) cos27b-u(pb) sin27B-s(pB) 

(17) 

and 

d2 = i(PA)-q(PA) cos27a-u(j9a) sin27A-s(pA)-i(j5B)+q(pB) cos27b+u(j9b) sin27B+s(j9B). 

(18) 

Note that for the combination of radiometers outputs corresponding to Stokes I, di + d2, the 
spurious term cancels, and for the combination corresponding to polarization signals, di— d2, 
it appears as a scalar map, independent of polarization angles 7a and 7b. Equation ([9]) is 
easily generalized to handle these extra terms by expanding the mapping matrix to 2A^t x4A^p, 
where the extra columns now correspond to the spurious signal map. Each row of this matrix 
now has 8 non-zero entries, 2 corresponding to each type of map, i, q, u and s. 



Input Transmission Imbalance 

An ideal differential radiometer only responds to differential input signals and completely 
rejects common mode signals. Due to several effects, the WMAP radiometers exhibit a small 
response to common mode signals (jJarosik et al.ll2003bl ). This response is characterized by a 



transmission imbalance factor, a^im, with the response of the radiometer, d, to input signals 
Ta and Tb given by equation [5l The values of Xim for the 3-year data are given in Table O 
The effect can easily be incorporated into the maximum likelihood maps solution (see eq. [9]) 
by modification of the mapping matrix, M. The TOD, including both the spurious map 
term and the transmission imbalance terms is given by 



di = (1 + Xini){i{pA) + q(PA) COS 27a + u{pa) sin 27a + s(pa)) + 
(1 - Xi^){-i{pB) - q(pB) cos 27b - u{pb) sin 27b - s(pb)) 



(19) 
(20) 



and 



'1 + Xim)(i(pA) - q(PA) COS 27a - u{pa) sin 27a - s{pa)) + 
(1 - Xim)(-i(j5B) + q(pB) COS 27b + u(pb) sin 27b + s(j9b)). 



(21) 
(22) 



The corresponding modification to M is that each term involving A-side data is mul- 
tiplied by a factor (1 + xim) and those involving B-side data are multiplied by a factor 



:i 
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3.4-2. Map Evaluation 

The maximum likelihood map solution (eq.[9]) is evaluated in two steps. First the 
product to = M"^N~^d, the "iteration 0" maps, are formed. The to maps are then multiplied 
by the (MN^^M^)^^ term using an iterative technique. 



Generation of the tg Maps 

The calibrated TOD, d, represents the entire sky signal, including the CMB dipole compo- 
nents. To minimize numerical errors a TOD signal corresponding to a nominal CMB dipole 
is removed from the calibrated TOD before evaluation of the to maps. The dipole signal 
removed from the TOD for each radiometer includes the effects of loss imbalance based on 
the values given in Table [H A CMB dipole amplitude of 3.3463 mK (thermodynamic) in 
the direction (/, b) = (263?87, 48?2) is used to calcul ate the barycentric c omponent. A CMB 



monopole temperature of 2.725 K (thermodynamic) (iMather et al.lll999l ) is used to calculate 
the annually varying component due to the spacecraft's motion about the barycenter. The 
kinematic quadrupole in not removed from the TOD in the generation of the to maps. 

The dipole-subtracted TOD from each radiometer must be multiplied by a filter con- 
sisting of the inverse of each radiometer's noise correlation matrix. This is accomplished by 
forming the inverse noise correlation function. 



lAtl < At„ 



lAtl > Atr 



where N{t') is the parametrized noise correlation function as described in § 12. 5[ Since N(t') 
is constructed to be zero at lags greater than Atmax, A^tt^ ^^^o falls to zero on the same time 
scale, so values of the filter, N^-^, are also set to zero for lags at which N{t') was set to zero. 
An example of the steps involved in forming this filter function for the Wll radiometer is 
shown in Figure [51 Next, a constant K is added to the values of N^-^ for lags less than 
^^max to force the mean of this portion of the filter to zero. This ensures that very low 
frequencies, those with periods longer than the filter extent, are filtered out. Finally, the 
filter is normalized by an overall multiplicative constant, C, such that the value at At = is 
unity. Note that none of these adjustments to N^'^ biases the resultant sky maps since the 
same form of N^'^ is used in both terms in equation (Q. Tabulated versions of these filters 
are contained in the data release. 

The radiometer noise properties are treated as stationary for each year of data, making 
the inverse noise matrix, N~^, circulant. The matrix multiplication, ■ d, may therefore 
be implemented through use of the convolution N^^ d, and is implemented using standard 
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Fourier techniques. Missing samples in the TOD, due either to gaps in the data or mask- 
ing when producing spm maps, are zero-padded to properly preserve the time relationship 
between data on either side of the gaps. 

Given the filtered TOD, the to maps are evaluated by accumulating the product of the 
corresponding elements of and the TOD sample by sample. Cut sky maps are formed 
by simply replacing the rows of M with zeros for observations where the data are masked. 
These tg maps correspond to the sky maps weighted by the inverse of the pixel-pixel noise 
correlation matrix, 

to = M'^N-M = M^N-^Mt = S^H (24) 

where t represents sky maps of the three Stokes parameters, I, Q and U, and a map cor- 
responding to the spurious signal described in § 13.4. 1[ These maps are available with the 
data release. Solving for the sky maps simply requires multiplication by the pixel-pixel noise 
covariance matrix, S. 



Final Sky Map Production 

The final step in the sky map production, multiplication of the to maps by the pixel-pixel 
noise correlation matrix, S, is effected through a conjugate gradient iterative technique. 
This method allows solution of the symmetric positive definite system 

Ax = b (25) 

for vector x given vector b and the ability to multiply an arbitrary vector by the matrix 
A. In this application vector b is a to map and matrix A corresponds to the inverse of the 
corresponding pixel-pixel noise correlation matrix, A = = M-'^N^^M. Multiplication 
of this matrix by a vector representing a set of sky maps is a straightforward operation. A 
simulated time stream is generated from the input maps by stepping through the archive 
evaluating the 16 non-zero elements of the mapping matrix (contained in two rows) corre- 
sponding to each observation (the M operation). The resultant time stream is then filtered 
and re- accumulated into a new map (the M^N~^) operation using the same method pre- 
viously described to produce the to maps. The conjugate gradient method then iteratively 
improves estimates of x, with the level of convergence being characterized by the quantity e 
given by the relation 

||b-Ax|| ||to-S-H|| 

e = [TTT, = ^ , (26) 

l|b|| ||to|| 

where is the vector-norm operator. The iterative processing was stopped when the 
condition e < 10^^ was satisfied. 
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Preconditioning 



The rate of convergence of the conjugate gradient method is quite sensitive to the properties 
of the matrix A, with nearly diagonal forms preferable. If another tractable matrix mul- 
tiplication operation can be found which approximates multiplication by A~^, the rate of 
convergence can be greatly increased. Consider the case where the matrix meets these 
criteria, multiplication of equation fl25l) by A^^ yields 



which has the same solution as equation fl25l) . but the matrix product A^^A is nearly 
diagonal, speeding convergence of the solution. This technique was used in production of 
the final sky maps. The preconditioner was implemented by separating the input map into 
high resolution (r9) and low resolution (r4) components. The low resolution map was formed 
by re-binning the high resolution map with uniform weights. This low resolution map was 
then subtracted from the high resolution map, leaving only small angular scale information 
in the high resolution map. The low resolution component was multiplied by the inverse 
of pixel-pixel noise correlation matrix (§ 13. 5p and the high resolution component multiplied 
by the reciprocal of the diagonal components of the full resolution noise correlation matrix. 
The two components were then recombined to produce the preconditioned map. The map 
production algorithms were carefully checked using numerous simulations to verify that they 
converged to the correct solutions. It was also verified that this solution was independent of 
the exact form of the preconditioner employed, although as expected the convergence rates 
did vary. Using the preconditioner as described, the r9 sky map sets, comprising Stokes I, Q 
and U components, and the spurious map S, converged in 50 — 100 iterations to the level 
described above. 



Both full sky and cut sky versions of maps were produced for all three years of obser- 
vations. The final sky maps for each year were produced by using data from the full sky 
maps to fill in the regions of missing data in the masked maps. Before combining the data 
from these maps the mean temperature of the full sky Stokes I and the S map required 
adjustment. For an ideal differential radiometer the mean value of the I and S maps would 
be undetermined, corresponding to singular modes in S. However, inclusion of the non-zero 
valued transmission imbalance factors, Xijn, converts these previously singular modes into 
modes with very small eigenvalues, indicating that they are very poorly constrained by the 
measurement. While no physical significance is attached to the recovered values, allowing the 




(27) 



3.4-3. Combining Full Sky and Cut Sky Maps 
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mean to vary is required to achieve the level of convergence previously described with regard 
to the conjugate gradient solution. The fact that the map means are poorly constrained 
indicates that estimates from the two different map processings may differ significantly due 
to statistical fluctuations. 

Failure to adjust the relative values of the means would introduce an obvious discon- 
tinuity at the border of the masked regions when the maps were combined. The procedure 
used to combine the maps is as follows. First the set of pixels which have the same number 
of observations in both versions of the sky maps were identified. A constant was then added 
to the full sky map such that the mean values of the previously identified sets of pixels in the 
full sky map equaled the mean value of the same set of pixels in each corresponding cut sky 
map. The pixels containing no observations in the cut sky map were then set to the values 
of the corresponding pixels in the adjusted full sky maps for the I and S maps. The mean 
values of the Q and U maps are well determined and were not adjusted. This procedure 
preserves the mean value for all four map components (I, Q, U and S) from the cut sky 
maps. Since the matrices were calculated using the cut sky coverage, cut sky maps 
with the appropriate mean value for use with these matrices may be obtained by simply 
masking the Galactic plane region. 

3.5. Evaluation of the Inverse Pixel-Pixel Noise Matrix, 

Several steps in the WMAP data processing and spectral processing require an explicit 
numerical representation of the inverse pixel-pixel noise correlation matrix, S^^. Formally 
this matrix is described as 

= M^N-^M = S-i(pi,p2) = Y,M{tuPi)N-,\U-t2)M{t2,P2), (28) 

where and M are the inverse noise matrix and the mapping matrix as described in 
§ 13.4.21 and § 13.4. 1^ and pi and p2 are pixel indices spanning the I, Q, U and S maps. The 
sums over ti and ^2 extend over all non-zero values of the inverse noise matrix evaluated at 
time ti — t2- These matrices were evaluated at r4 on a year- by year basis for all 10 DAs, 
resulting in 30 matrices. One approximation was used in the evaluation of these matrices to 
speed processing. Recall that each observation populates two rows of the mapping matrix, 
each row corresponding the projection factors for each of the two radiometers comprising 
the DA. Within each row the 8 non-zero elements correspond to the two differential pixels 
(A-side and B-side) of the 4 maps I, Q, U and S. The approximation consists of grouping 
time contiguous observations for which both the A-side beam and B-side beam remain in 
the same r4 pixels. Each such group comprises ~ 30 rows of the mapping matrix (about 15 
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observations) whose non-zero elements reside in the same 8 columns. For each contiguous 
group of observations, the averages of each of the columns containing non-zero values is 
calculated for the rows associated with each radiometer. These averages are tagged with 
the group's starting and ending times, the pixel indices of the A-side and B-side beams 
and the radiometer identification. Given the starting and stopping times it is possible to 
calculate the appropriately weighted sum of N^^{ti — ^2) corresponding to any pair of groups 
and corresponding radiometer. The sums over ti and t2 are then calculated for all pairs of 
groups for which the weighted sum of N^^{ti — t-i) is non-zero. This is accomplished by 
multiplying this weighted sum by the average values of the appropriate columns of the 
corresponding groups, and accumulating the values into pixel-pixel noise matrix elements 
based on the pixel indices associated with each group. The result is a 4A^p x 4A^p matrix, 
S^^, which may be inverted by standard numerical techniques to form the noise matrix S. 
Using the group average values for the data in the columns of M rather than the exact row 
by row values in calculation of the inverse noise matrix is a good approximation for the vast 
majority of the sky, since the values being averaged are smooth and slowly varying. However 
these assumptions fail for pixels very close to the Galactic poles since the polarizations 
angles 7a and 7b vary significantly within a r4 pixel. The effect of these approximations 
has been investigated and is found to have negligible affect on the calculated power spectra. 
The r4 noise matrices contained in the data release are calculated from mapping matrices 
corresponding to cut sky coverage and are intended for use in regions outside of the region 
excluded by the processing mask. Values corresponding to r4 pixels entirely contained within 
the processing mask are set to zero. 



3.5.1. Projecting Transmission Imbalance Modes from the S ^ Matrices 

One potential source of systematic artifacts in the sky maps arises from the error in the 
determination of the input transmission imbalance parameters, ,Xim. These parameters are 
measured from the fiight data and arc used to calculate the estimated dipole signal that is 
removed from the TOD in preparation for the conjugate gradient map processing. Due to the 
large amplitude of the dipole signal, even small errors in the measured values of transmission 
imbalance parameters could introduce significant artifacts into the sky maps. This effect was 
studied through the use of simulations in which a simulated TOD archive was produced with 
a given set of values and was analyzed with the input value of as well as values 20% 
higher and lower to simulate errors in the measurement of x^^. Differences between the 
maps produced using the correct and incorrect displayed well defined spatial structure 
confined to low I with small variations between DAs arising from the slightly different scan 
patterns. The geometry of these structures is completely determined by the scan patterns. 
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while their amplitude in the final sky maps is determined by the difference between the 
value used to process the TOD and the "true" value of Xim- The patterns from these 
simulation are therefore treated as templates, used to exclude the modes corresponding 
to the aforementioned spatial structures from subsequent analysis. This is accomplished 
through the use of modified versions of the inverse noise matrices, S~^, that have these 
modes projected out. Each modified matrix is calculated as 



so that 



0. 



(29) 



(30) 



In the above expressions designates an outer-product and v is the map mode to be removed, 
obtained from differences between the simulated maps made with the correct transmission 
imbalanc e factors and the incorrect factors. Using these forms of the matrices in the spectral 
analysis (jPage et al.l l2006l ) eliminates artifacts which might result from the errant signals 
cause by errors in the measured values of the transmission imbalance parameters. Figure [T2 
shows the effect of suppressing these modes on the low-^ polarization power spectra. This 
correction has a very sm all effect on the c osmological parameters, but has been included in 
the results presented in iPage et al.l (120061 ) for completeness. 



Inverse noise covariance matrices based on these correlation matrices are contained in 
the data release. The covariance matrices are the correlation matrices scaled by a multi- 
plicative factor, I/ctq, where ctq is the noise per observation corresponding to each DA. 



3.6. Production of Low Resolution Maps 

Evaluation of the low-i power spectra with optimal signal-to-noise involves use of sky 
maps and the inverse pixel-pixel noise matrices, S~^, which describe the maps' noise proper- 
ties. Since the matrices are produced at r4 due to computational constraints, corresponding 
low resolution sky maps must be produced that preserve the signal and noise properties as 
well as possible. There are several difficulties involved with production of low resolution 
maps. Simply binning maps to lower resolution, using either uniform or pixel-by-pixel (di- 
agonal) inverse noise weighting will result in aliasing of high-£ signal and noise into lower-£ 
modes. This effect can be reduced by applying a low-pass filter to the map before degrada- 
tion, but such a filter introduces additional noise correlations, not described by the inverse 
noise matrix as described in the previous section. Even without application of a filter, 
the omission of the inter-pixel noise correlations results in the degraded maps having noise 
properties not precisely described by the noise matrix as calculated in § 13.51 
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Fig. 12. — Difference between the low-i polarization power spectra witli and witliout removal 
of the transmission-imbalance modes from the matrices. The black and red lines are the 
differences between the power spectra of combined Q-band and V-band polarization data for 
the E and B modes respectively. For comparison, the dashed line is the expected E-mode 
signal corresponding to the best fit WMAP power spectrum. The effect is largest for the 
B-mode i = 3 which already has a large uncertainty arising from the scan pattern 
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It is possible to produce low resolution sky maps directly using the procedures described 
in § 13.4.21 by accumulating data directly into low resolution maps. This corresponds to 
reducing the number of columns in the mapping matrix, M, to reflect the reduced number 
of pixels in the maps. Such maps will have noise properties accurately described by the 
pixel-pixel noise matrix, but will have a higher degree of signal aliasing than maps initially 
produced at high resolution and subsequently degraded. This occurs since in evaluation of 
the sky map (see eq.[9]) the low resolution form of the mapping matrix is apphed to the data 
multiple times, each time introducing aliased signal components, whereas degrading a map 
initially produced at high resolution only introduces the aliased components once. 

For the Stokes I maps, which are highly signal dominated, degraded res 9 maps are 
clearly preferable, since the small errors in the noise description are of no consequence and 
these maps contain lower levels of aliased signal. For the Stokes Q and U maps the signal to 
noise ratio in the low-£ regime is on order unity. Low resolution maps of Stokes Q and U were 
produced both directly at r4 and by degrading r9 maps without application of an anti-alias 
filter. Analysis of both sets of maps produced very similar low-i polarization results for both 
simulated and flight data. Based on these results it was determined that the noise properties 
of the degraded form of low resolution polarization maps was adequately described by the 
pixel noise matrices, so the degraded form is used for all three Stokes parameters. 

Although the inter-pixel noise correlations are ignored in the map degradation, the 
noise correlations between the Stokes Q and U signals within each high resolution pixel are 
included. Two different versions of reduced resolution maps are produced and are described 
in the following sections. 

Production of Full Sky Reduced Resolution Maps 

The full sky reduced resolution maps (r4) are comprised of two components. The first 
component, which excludes the Galactic plane region, is obtained by degrading the r9 cut 
sky (spm) maps and has noise properties described by the inverse noise matrices ( §3.5p . The 
second component, used in the Galactic plane region, results from the degradation of the 
combined maps ( §3.4.3p and has noise that is not described by the inverse noise matrices. The 
noise in this portion of the maps is approximately described by the number of observations 
values, Nohs, contained in the sky map data products. 

The first component of the reduced resolution map, comprising the high Galactic latitude 
regions, is generated as follows: During production of the high resolution spm maps, the 
total weights applied to each pixel in the I, Q, U and S maps, iVobs, Nqq, Nuu, and N^s, 
are calculated. These are simply the sums of the squares of the corresponding columns of 
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the mapping matrix, M. The QU, QS and SU off-diagonal terms, A^'qu, A^'qs and A^'su, with 
each r9 pixel are also calculated. These are the sums of the products of the corresponding 
columns of the mapping matrix. Figure US] shows maps of the Nohs, Nqq, Nuu, and Nqu 
patterns. 

For each high resolution pixel, p9, an Nobs matrix is formed. 





f ^QQ 


Nqv 
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Nobs(p9) = 1 
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(31) 
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which differs only by a multiplicative constant from the inverse noise matrix. Correlations 
between the Stokes I and the polarization related components are not used in the map 
degradation. The Q, U and S values in a low resolution map pixel, p4, are then simply the 
inverse noise weighted average of all the Q, U and S values of the high resolution pixels it 
contains, 

(Qp4 \ I Qp9 \ 

UpA = (N*,°t.s(p4))"' E ^(p9)Nobs(p9) ^p9 , (32) 

Nrbs(p4) = J] «;(p9) N„,,(p9). (33) 

p9Gp4 



Here w{p9) is the r9 processing mask (IHinshaw et al.l l2006l ) containing values of and 1 



and the summations are for all r9 pixels contained within each r4 pixel. Elements of the r4 
maps containing one or more observations are retained. The Stokes I maps were degraded 
using simple Nohs weights. 

Elements of the r4 maps described above with no observations are filled with values from 
degraded forms of the combined maps ( §3.4.31) . The combined maps are degraded following 
the same procedure, but the Nobs(p9) matrices are formed using weight values from the cut 
sky (spm) maps when outside the processing mask, and values from the full sky (fs) maps 
within the processing mask. Since the weights corresponding to the full sky maps differ from 
those of the spm maps, the inverse noise matrices, calculated from spm map weights, do not 
describe the noise in this portion of the map. All the rows and columns of the matrices 
corresponding to these pixels contain zeros. 



Production of the Partial Sky, Foreground Cleaned Reduced Resolution Polarization Maps 



The r educed resolution sky maps used in the pixel based low-£ likehhood evaluation (|Page et al 



20061 ) are formed from the foreground cleaned Q and U r9 spm maps, as described in lHinshaw et al 
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Fig. 13. — Maps of Nohs, Nqq, Nuu, and Nqu weights for K band year-1 spm sky coverage. 
These maps are measures of the noise covariance and are used when the maps are degraded 
from r9 to r4. 
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( 120061 ) and iPage et al.l (120061 ) . Including the spurious term in the map degradati on and hke- 
lihood evaluation was found to have a negligible effect on the likelihood results (jPage et al. 
20061 ). Therefore, to simplify and speed calculation, only the Q and U terms were used in the 
the d egradation and low-i likelihood evaluation of the foreground cleaned maps (jPage et al. 
2006h . 



The foreground cleaned Q and U maps were degraded in a similar fashion to the full 
sky maps described above: For each high resolution pixel, p9, an Nobs niatrix is formed. 



(34) 



The values of the low resolution map pixels, p4, are the inverse noise weighted sum of 
the values of all the constituent high resolution pixels. 



= mp4)(Nrs(p4))-^ 5:«;(p9)No,,(p9)(^P^ 

^ p9gp4 ^ 



Nrs(p4) = 5^ ^(p9)N„,,(p9). 

p9Gp4 



(35) 
(36) 



In this case w{p9) is the r9 P06 mask (jPage et al.l 120061 ) containing values of and 1 and 
Vl^(p4) is a low resolution form of this mask which contains zeros in elements for which more 
than half of the r9 pixels were masked and values of unity otherwise. These degraded maps 
are contained in the data release. 



4. TESTS ON THE DATA SET 

4.1. Limits on Spin-Synchronous Artifacts in the TOD 

Among the most troublesome type of systematic errors in experiments such as WMAP are 
spurious signals which are synchronous with the scanning motion of the instrument. Such 
artifacts do not integrate to lower levels with the inclusion of additional data as does random 
noise. The most likely source of such driving signals is from varying insolation due to the 
observatory's motion driving temperature variations in the instrument hardware, voltage 
fluctuation on the spacecraft power bus, or coupling directly into the microwave optics. For 
all of these effects the errant signal is expected to be synchronous with the 129 s motion 
of t he sun relative to th e spacecraft. Following the procedure used in the year-1 analy- 



sis (jjarosik et al.l l2003bl ) the radiometer output signals were binned in a spacecraft fixed 
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coordinate system based on the azimuth of the Sun about the spin axis of the spacecraft. 
When binned in such a manner any signal synchronous with the azimuthal position of the 
sun about WMAP should be preserved while asynchronous signal components should average 
away. Combinations of radiometer outputs corresponding to both temperature anisotropies, 
dl + d2, and polarization signals, dl — d2, were accumulated for the entire 3 year mission. 
The resultant signals were all found to be less that 1.25 /xK rms. We therefore believe that 
there is little contamination of the TOD arising from spin synchronous solar effects and no 
corrections are applied to the TOD for spin synchronous signals. 

4.2. Year-to-Year Null Tests 

One of the most fundamental tests of the WMAP maps is the year-to-year consistency 
between maps. Although such tests do not eliminate the possibility of processing or system- 
atic errors that repeat year to year, they do test for spurious signals such as striping from 
radiometer 1// noise and glitches. A comparison between the year-1 and 3-year processing 
of the first year data from the V2 DA was presented in § 12. 4. H and displayed a predominantly 
quadrupolar feature associated with the use of the improved radiometer gain model. Since 
both maps used in that comparison were generated from the same TOD, the white noise 
component largely canceled in the difference maps allowing this feature to be seen. In the 
case of year-to-year comparisons such white noise cancellations do not occur, so subtle effects 
may not be evident. Nevertheless two sets of year-by-year difference maps are presented to 
demonstrate the year-to-year consistency, one based on K bands with the largest foreground 
signal, and one based on V-band with the smallest foreground signal. The top four panels 
in Figure [H] display the 3 year average and year by year differences of the K-band full sky 
maps. The 3-year combined map is dominated by Galactic emission and also displays a num- 
ber of point sources at high Galactic latitudes. The high level of Galactic emission makes 
K-band one of the more difficult bands to calibrate. While the year-by-year null maps are 
noise dominated some residual Galactic plane signal is evident. This signal arises from small 
calibration errors of ~ 0.1 — 0.2% and is consistent with the stated calibration uncertainty 
of 0.5%. Also evident in these maps are some features of both signs at high latitude arising 
from known variable sources. 

The lower four panels in Figure [T3] display the 3- year average and year- by-year difference 
of the V-band data. In this case Galactic emission can be seen at low latitudes and CMB 
anisotropies at high latitude. The difference maps are again noise dominated with the 
spatially varying noise properties evident. The increased noise near the ecliptic plane is a 
consequence of the lower number of observations in these regions as determined by the scan 
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strategy. 

Simple power spectra of individual year-to-year null maps are dominated by the ra- 
diometer noise bias which grows as ^{l + 1) when plotted with the customary scaling, and 
are not useful for assessing year-to-year map consistency. Figure [15] compa res the best fit 



WMA P CMB power spectrum with the year-1 to year-2 cross power spectrum (IHinshaw et al. 



2003bj) of combined V and W band Stokes I sky maps, averaged in A£ bins of 100. For each 
year, the linear combination of sky maps (V1-V2) -|- (W1-W2) -|- (W3-W4) was formed. 
This combination removes the true sky signals to the extent that the bandpasses and gains 
of each pair of maps are matched, while the use of the year-to-year cross spectrum eliminates 
the noise bias since the radiometer noise is uncorrelated year-to-year. The power spectrum 
of the null maps are very small relative to the measured power spectrum. 

A set of null maps for the Stokes Q and U parameters is displayed in Figure [12] for 
K-band. Here in the year-combined maps the Galactic foreground emission, predominantly 
synchrotron, is clearly evident. The year-by-year difference maps are again noise dominated. 
The spatially varying noise structure for the Stokes Q and U maps is more complicated than 
that of the Stokes I maps due to the additional polarization projection factors (e.g. sin 27) 
and the noise covariance between the Stoke Q and Stokes U maps of each DA. The faint cross- 
like features near the ecliptic poles are regions with large weights, and therefore low statistical 
noise, as can be seen by comparing these null maps with the iVobs maps in Figure [13] Simple 
null tests, as displayed in Figure [T^] for Stokes I, are not useful in assessing year-to-year 
consistency of the polarization maps due to the much smaller level of the expected signals 
and more complicated noise properties of the polarizat ion maps. Detail ed analysis of the 



polarization results, including null tests, is presented in [Page et al.l (l2006l ). 



4.3. Determination of Map Noise Levels 

Evaluation of the instrument noise level is an important aspect of many analyses involv- 
ing the WMAP sky maps. In principle, a single noise per observation parameter, ctq, should 
be sufficient to characterize the noise level for each set of maps (a set consisting of the I, Q 
and U maps for a given year of data and DA). Such a characterization would require the 
use of full resolution pixel-pixel noise correlation matrices, describing the noise correlations 
between all combinations of pixel indices and Stokes components. For r9 sky maps, the size 
of such matrices precludes their use. However, it is possible to characterize the map noise 
levels using only the diagonal components of the noise correlation matrices, A^obs, A^qq, and 
Auu, provided the intra-pixel components of noise correlation matrix, A'qu, A'qs, and A'su 
(§ 13.61) are included. This approximation uses different values of the noise per observation 
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K1 3-year Average K1 year1 - year2 




Fig. 14. — WMAP 3-year combined and year-by-year difference maps of Stokes I for K- 
band and V-band. Tlie residual Galactic plane features in the K-band difference maps are 
consistent with the absolute calibration error of 0.5 %. Several variable sources are also 
visible in the K-band maps. The CMB anisotropy signal is visible in the combined V- 
band maps while no signal is visible in the V-band year-by-year difference maps. Maps are 
displayed at r5. 
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Fig. 15. — Cross power spectrum of year-1 and year-2 null maps. For each year Stokes 
I null maps were formed for DA combinations ((V1-V2) + (W1-W2) + (W3-W4))/6, and 
the year-to-year cross spectrum calculated. Within each year the map differences remove 
true sky signals, and the year-to-year cross spectrum eliminates the noise bias. The power 
spectrum (diamonds) has been binned in A£ — 100 to reduce noise. For comparison the solid 
line is the best fit WMAP Stokes I power spectrum. The null map residuals are very small 
relative to the measured power spectrum, reflecting the repeatability of the measurements 
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Stokes Q 3-year Average Stokes Q year1 - year2 




Fig. 16. — WMAP 3-year combined and year-by-year difference maps of Stokes Q (top 4 
panels) and U (bottom 4 panels) for K-band. Galactic foreground emission is clearly visible 
in both Stokes Q and U of the three year combined maps. The difference maps show no 
residual signal, although the spatially varying noise levels are apparent. Maps are displayed 
at r5. 
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parameter for maps of the different Stokes components, parameterizing the effects of the 
ignored, off diagonal components of the noise correlation matrices. The different values of 
these parameters for the different map components result from the differences between the 
frequency distributions of the intensity and polarization information in the TOD, and the 
radiometers' frequency dependent, i.e. 1/f, noise properties. The noise levels of the radiome- 
ters are also found to be stable year-to-year, so a single noise per observation value for each 
sky map/DA combination describes all three individual year sky maps. 



4.3.1. Noise Evaluation of the r9 Maps 

The noise levels of the r9 maps in the current data release have been estimated using 
two different methods. In the first method the noisc-pcr-obscrvation, is calculated directly 
for each DA from a weighted average of the variance of year to year differences of sky maps. 
For the I maps the noise per observations based on year x and year y sky maps, ctq^'^ is 
given as 



A^pix = (38) 
p 

where t'^{p) and n^-^^{p) are the Stokes I value and the number of observations for pixel p 

of the year x sky maps, and w(p) is the processing mask, used to exclude regions of high 
Galactic emission from the analysis. Extending this technique to the polarization maps 
requires a generalization to allow for the noise covariancc between the Stokes Q and U maps 
and the spurious map, S. This is implemented using the relations 

= N^^Y.''{p).{z^y{p)r {M^y{p)r' z-y{p) (39) 
p 

(Ql-Ql\ 

z-y{p) = U--UI (40) 

J^'^ip) = NS,,(p)-i + Nl(p)-^ (41) 

where Z^y(p) is a vector comprised of the differences between the year x and year y values 
of the Q, U and S maps at pixel p, and N^jjg(p) is the corresponding 3x3 element matrix 
describing the intra-pixel noise correlations between the Q, U and S maps for pixel p and 
year x. Relative calibration errors between the two maps and the year to year differences 
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in the distributions of observation pointings within pixels will allow some sky signal to leak 
into the difference maps, biasing the resultant noise estimates high. Values of ctq and ctq 
obtained by averaging the values of the three possible year-by-year difference combinations 
are presented in Table [5l 

The second technique also measures the noise from the year to year difference of the sky 
maps. Difference maps of each Stokes parameter, I, Q and U, are formed from year x and y 
sky maps for each DA. The effective number of observations for each pixel of these difference 
maps are calculated. For the Stokes I the effective number of observation is given as 

n:^\p) = 1/ {KMr' + KM)-') (42) 
and for the Stokes Q and U components it is given as 

<^^(p) = (43) 

and 

<^'"(p) = 1/Klip) (44) 

where and J\f^^ are the diagonal components of M^^ip) as defined in equation (14T1) . For 
each Stokes component the following procedure is applied: 1) Pixels are divided into groups, 
Gi, based on their values of riefr, 2) The variance of the difference maps for the pixels in each 
group, Var(G'i), are calculated as well as the mean value of iVeg, rieaiGj), 3) These values are 
fit with two parameters, and a^^ in the form 

Var(GO = {al^ f ■ ^^:^) + {a^^f (45) 

where a^^ is a parameter that allows for the presence of a statistically isotropic signal that 
does not scale with neff(Gi) . Pixels in high Galactic emission regions, as determined by the 
processing mask, are excluded from the analysis. Fits of all three year-by-year difference 
maps are performed and the resultant values averaged together, these results are also pre- 
sented in Tableland agree within 0.3% with the values obtained using the first technique. 
Since there is no significant disagreement between the methods we use the results of the first 
technique for analysis. By comparing the noise level obtained from the year 3 - year 2 dif- 
ference maps with those from the year 2 - year 1 difference maps we find that the sensitivity 
of the radiometers has decreased at most by 0.1% over the 1 year period between the mean 
epoch of the data sets. 

4-3.2. Noise Evaluation the r4 Full Sky Maps 

The noise of the full sky reduced resolution (r4) maps may be determined using the 
pixel-pixel noise matrices. This analysis is performed only on the high Galactic latitude 



- 53 - 



Table 5. Noise Levels for r9 Maps 





Method 


1 


2 


1 


2 


2 




DA 


Signal 


-I 












Kl 




1.439 


1.439 


1.455 


1.453 


1.453 


1.057 


Kal 




1.464 


1.462 


1.483 


1.484 


1.480 


1.019 


Ql 




2.245 


2.244 


2.269 


2.273 


2.272 


1.027 


Q2 




2.135 


2.140 


2.156 


2.157 


2.156 


1.016 


VI 




3.304 


3.311 


3.330 


3.333 


3.333 


1.023 


V2 




2.946 


2.951 


2.970 


2.970 


2.974 


1.015 


Wl 




5.883 


5.896 


5.918 


5.935 


5.930 


1.023 


W2 




6.532 


6.550 


6.571 


6.587 


6.578 


1.024 


W3 




6.885 


6.885 


6.925 


6.918 


6.931 


1.020 


W4 




6.744 


6.765 


6.780 


6.793 


6.800 


1.033 



Note. — The values presented in the first 5 columns of this table are estimates of the noise 
levels of the Stokes I, Q and U parameters in the 3-year high resolution maps for each DA. 
Each value was calculated by two different methods as described in § 14.31 and the results are 

1 /2 

presenting in units of mK ■ A''^^ • The two methods agree to within 0.3% for both the high 
resolution Stokes I and polarization maps. The values in the last column are the ratios of 
the low resolution (r4) maps noise estimates, obtained using the noise correlation matrices 
( §4.3.2p . to the high resolution noise of the I maps ( §4.3.ip . Differences in the noise estimates 
arise from ignoring pixel-pixel noise correlations both in the map degradation and the noise 
estimation of the high resolution maps. 
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regions, those pixels originating from degradation of the high resolution spm maps. For each 
DA, noise levels are evaluated for all three possible year-to-year difference map combinations. 
For combinations of years x and y the noise per observation, ctq'^ , is calculated as 

(t.-t,)^S-J(t.-ty), (46) 
(S. + Sy)-\ (47) 

where t^; and ty are reduced resolution sky map sets comprising I, Q, U and S components for 
years x and y respectively, and and Sy are the pixel-pixel noise correlation matrices for 
the corresponding years. For each DA, the values from the three year-to-year combinations 
are averaged together. The ratio of these average values to the noise values of the r9 I maps 
are presented in Table [51 The differences between the two noise estimates arise from ignoring 
the pixel-pixel noise correlations in both the r9 - r4 map degradation, and the evaluation of 
the (To of the r9 maps. 

5. SUMMARY 

WMAP continued to operate normally throughout its second and third years of science 
observations. The additional observations served not only to reduce the statistical noise, 
but also provided for improved understanding of the beam geometries and radiometer per- 
formance, important factors in interpretation and processing of the science data. Updated 
versions of the radiometer gain model, beam models and window functions were presented. 
New processing procedures that produce maximum likelihood sky maps for Stokes I, Q and 
U parameters have been described in detail as well as the calculation of the inverse noise 
matrix. The techniques used to quantify systematic errors were described along with details 
of the sky map processing used to correct for potential systematic artifacts. Sky maps pro- 
duced from data from different epochs were used to display the year to year reproducibility 
of the signal and to measure the instrument noise levels. The performance of the WMAP in- 
strument remains nominal with no significant decrease in sensitivity observed to date. 

The WMAP mission is a collaborative effort between Goddard Space Flight Center, 
Princeton University, the National Radio Astronomy Observatory and is supported by the 
Office of Space Sciences at NASA Headquarters. EK acknowledges support from the Alfred 
P. Sloan Fellowship. Some of the results in this paper have been derived using the HEALPix 
(Gorski, Hivon, and Wandelt 1999) package. 
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A. Sidelobe Correction and Rescaling 



Following a notation similar to that of iHinshaw et al.l (l2003a] ) the uncalibrated time 



ordered differential data for each detector, c, may be represented as 

c = sj rffi„M(n)t(n). (Al) 

Here g is instrument responsivity (cts/mK), M(n) is the mapping function which transforms 
spatial information into temporal information and t(n) is the sky signal being observed (mk). 
All the aforementioned quantities are time dependent as indicated by the use of boldface 
symbols. 

The mapping function can be expressed in terms of the differential beam response of the 
instrument in instrument coordinates, i?d(ninst), and R, the time dependent rotation matrix 
that transforms sky coordinates, n, into the satellite coordinates, ninst- 

M(n) = 5d(R ■ n) (A2) 

The differential beam response is a signed quantity and is normalized such that 

rff]n|5d(n)| = 2. (A3) 



4tt 

The fraction of the differential beam response in the sidelobes, /si, is obtained by integrating 
the magnitude of the differential beam response over the region of the map further than th 
from the boresights of both main beams, 

/si = i/ dnr.\Bf{n)l (A4) 



2 



dA,rfB>m 



l eaving the fraction (1 — /si) of the gain in the main beams. As described in lHinshaw et al. 



(j2003al ). the calibration of WMAP is tied to the amplitude of the annual dipole signal arising 
from the motion of WMAP with respect to the solar system barycenter. If the annual dipole 
sky signal is represented as tad, the corresponding component of the uncalibrated radiometer 
signal is 

Cad = g J rf^]„M(n)tad(n). (A5) 

Calibration is accomplished by performing a least-squares fit of the uncalibrated data to a 
predicted time ordered data set, d!J^°°°, based on the known monopole temperature. 

Cad = g d~, (A6) 
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yielding an estimate of the radiometer gain, g. The predicted time ordered data set is 
obtained using a mapping function assuming ideal 'pencil' beams, MP'^"'^'' and the annual 
dipole sky signal, 

d™ = J rf(]nMP"""^V)tar°(n). (A7) 

The actual calibration procedure is iterative and provides for subtraction of better estimates 
of the fixed sky signal before each new gain determination. Once the iterative steps are 
completed the radiometer gain model (Sec. 12.41) is fit to smooth and interpolate the hourly g 
values. Taken together, the above cahbration procedure effectively measures g by evaluating 
the mean of the ratio 

g/rfr]„M(n)W(n) \ 
Jrffi„Mpencii(n)t™(n)/' 

Note that the mapping function in the numerator incorporates the full sky beam response. 
This function can be written as the sum of a main beam mapping function, M™*^, and a 
sidelobe mapping function, M*^\ 

M = M'"^ + M'K (A9) 
In terms of these quantities the estimate of the radiometer gain becomes 

- ^ / g / dQnM^\n)t,a{n) + g / dQnM^\n)t^a{n) \ 
^ \ /(il]„Mpencii(n)t™(n) /■ 



This can be approximated as 



ymono 

^1 /si "I" f /si) yjjjQj^Q • (-^H) 



The response of the telescope main beam to a dipole signal is nearly identical to an ideal 
pencil beam except that only the fraction (1 — /si) of the antenna gain is contained in the 
main beam, leading to the first two terms in the parenthesis in equation lAllI The last 
parenthesized term originates from the annual dipole signal picked up through the sidelobe 
response, with e being a coupling efficiency factor parameterizing how well the dipole signal 
couples to the sidelobe response pattern. The product e/si was evaluated by generating 
one year of simulated TOD representing the sidelobe response to a predicted annual dipole 
signal. This TOD was fitting in one hour periods to a template function corresponding 
to the expected pencil beam response signal. Although the quantity e is in principle time 
dependent, it was found to vary by less than 10% over the entire year simulation, and it is 
approximated as constant. Applying this calibration factor to the uncalibrated data yields 
a calibrated archive, 

d = ^ = [ ^^^]„M(n)t(n). (A12) 

g J- — /si + e/si J 
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A sidelobe corrected archive, calibrated to correspond to the main beam containing all 
the antenna gain can be written as 

d'"" = [ t/fi„M'"b(n)t(n), (A13) 

i- - J si J 

where the prefactor of the integral compensates for the fact that only the fraction (1 — /gi) 
of the beam response is contained in the main beam in M™^(n). To first order in f^i this 
can be obtained from the uncorrected calibrated archive, 

^ (1 + _ J dr]„M^V)t(n), (A14) 

where t(n) = tgxcd + is an approximate sky map obtained from applying the iterative 
map making algorithm to the calibrated (but not sidelobe corrected) TOD archive given by 
equation IA12I Values of the recalibration factor, (1 + efsi), are presented in Table IH 
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